# Street et al SSq 2015 replication.R

# created using:
# R version 3.1.0 (2014-04-10) -- "Spring Dance"

# set working directory, load data
# setwd("/Users/... YOUR DIRECTORY HERE .../")
load("streetetal_ssqdata_replic.RData")
ls()
# "a" is a data frame with the first set of imputed data, the main source of results
# "b", "c", "d" and "e" are four further sets of imputed data
# to check results in the other sets of imputed data...
# replace "a$..." below with "b$..." 
# replace "data=a" with "data=b"
# etc.
library(survey)
library(foreign)

###
### Removal statistics
###

# Starting in the year 2000, in thousands
removes <- c(188, 189, 165, 211, 241, 246, 281, 319, 360, 392, 382, 387, 418, 438)
# Sources: PewHispanic, US Dept of Homeland Security
# http://www.pewresearch.org/fact-tank/2014/10/02/u-s-deportations-of-immigrants-reach-record-high-in-2013/
# (last accessed Jan 6 2015)
mean(removes[2:9])
#[1] 251.5   so mean of 252k per year under Bush
mean(removes[10:14])
#[1] 403.4   so mean of 403k per year under Obama
403.4/251.5
#[1] 1.603976   so Obama deported at a rate around 1.6 times higher
# Survey wording: "In fact, the Obama administration has deported about one and a half times more people each year than the average under President Bush."



###
### Gallup presidential approval data (for figure 1 in the paper)
### 

### Source: http://www.gallup.com/poll/124922/presidential-approval-center.aspx

appro <- read.csv("HispObamApprov.csv")
# (note, this has most recent first, so need to reverse it)

want <- rev(appro$hispanics)
want2 <- rev(appro$blacks)
want3 <- rev(appro$whites)
length(want)

# records were available for download up to week of mar 17-23, 2014
# thereafter, no longer possible to download from gallup
# so add in the rest by hand

# (this takes me up to week ending april 26, 2015)

whisp <- c(want, 56, 57, 56, 59, 55, 60, 57, 54, 52, 51, 53, 50, 53, 55, 53, 52, 55, 50, 54, 52, 53, 55, 48, 44, 46, 48, 45, 57, 44, 49, 51, 55, 49, 54, 58, 68, 60, 68, 70, 61, 66, 63, 69, 65, 66, 62, 65, 65, 63, 67, 64, 61, 64, 64, 68, 66, 57)
length(whisp)

wblack <- c(want2, 79, 90, 82, 84, 84, 84, 88, 88, 86, 87, 87, 82, 82, 81, 86, 84, 87, 85, 85, 85, 86, 80, 89, 81, 83, 88, 85, 85, 81, 84, 87, 82, 84, 78, 86, 76, 82, 86, 84, 89, 92, 82, 83, 86, 85, 91, 92, 82, 85, 86, 88, 87, 89, 88, 91, 86, 86)
length(wblack)

wwhite <- c(want3, 33, 32, 33, 33, 34, 33, 34, 34, 34, 35, 34, 32, 30, 32, 33, 32, 32, 30, 30, 32, 30, 33, 30, 31, 31, 34, 33, 32, 32, 31, 31, 31, 29, 32, 31, 31, 32, 30, 32, 31, 34, 34, 33, 37, 36, 33, 35, 34, 33, 34, 35, 33, 34, 35, 35, 33, 34)
length(wwhite)


postscript(file = "figure1.eps", horizontal=F, onefile=F, paper="special", width=6.5, height=4.5, pointsize=9)

par(mar=c(2.1, 3.3, 2.3, 0.6))
par(mgp=c(2, 0.75, 0))
plot(whisp, ann=F, axes=F, ylim=c(27.5, 101), type="n")
axis(1, at=c(0, 49, 101, 153, 205, 257, 311), labels=c("2009", "2010", "2011", "2012", "2013", "2014", "2015"))
axis(2, at=c(20, 30, 40, 50, 60, 70, 80, 90, 100))
abline(h=20, lwd=0.5, col="grey")
abline(h=30, lwd=0.5, col="grey")
abline(h=40, lwd=0.5, col="grey")
abline(h=50, lwd=0.5, col="grey")
abline(h=60, lwd=0.5, col="grey")
abline(h=70, lwd=0.5, col="grey")
abline(h=80, lwd=0.5, col="grey")
abline(h=90, lwd=0.5, col="grey")
abline(h=100, lwd=0.5, col="grey")
lines(lowess(whisp, f=2/9), lwd=3, col="darkgrey")
points(whisp, pch=3)
lines(lowess(wblack, f=2/9), lwd=3, col="darkgrey")
points(wblack)
lines(lowess(wwhite, f=2/9), lwd=3, col="darkgrey")
points(wwhite, pch=2)
title(ylab="Weekly mean % of group approving of Obama's work as President", cex.lab=1)
legend(226.5, 104.5, legend=c("Blacks", "Latinos", "Whites"), pch=c(1, 3, 2), horiz=T, bg="white")
title(main="Obama job approval by race/ethnicity of respondent (source: Gallup)")
box()
dev.off()




############################################################################################
# Now the survey data: Code the main outcome variable
############################################################################################

table(a$demwelcom)
#         dk     neither     refused unwelcoming   welcoming 
#        228         281          15         119         407 
# (no imputation on this variable)

# prefer a version that treats "dk" as missing, and puts "neither" together with "unwelcoming"
# so the question is, welcoming or not?
a$demdv3 <- numeric(nrow(a))
a$demdv3[a$demwelcom=="welcoming"] <- 1
a$demdv3[a$demwelcom=="refused"] <- NA
a$demdv3[a$demwelcom=="dk"] <- NA
table(a$demdv3)
#  0   1 
#400 407

# so THIS is the core result
lm <- lm(demdv3 ~ toldtruth, data=a); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.55392    0.02466   22.46  < 2e-16 ***
#toldtruth   -0.10029    0.03507   -2.86  0.00435 **


## confirm that we get similar results in KN panel vs the online-recsruited

lm <- lm(demdv3 ~ toldtruth, data=a[a$DOV_PANEL=="KN PANEL MEMBER",]); summary(lm)
#(Intercept)  0.52632    0.04662  11.290   <2e-16 ***
#toldtruth   -0.11362    0.06434  -1.766   0.0787 .  
xtabs(~ a$toldtruth[a$DOV_PANEL=="KN PANEL MEMBER"] + a$demdv3[a$DOV_PANEL=="KN PANEL MEMBER"])
# so N is 114 ppl not told truth, and 126 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$DOV_PANEL=="OFF PANEL SAMPLE",]); summary(lm)
#(Intercept)  0.56463    0.02906  19.427   <2e-16 ***
#toldtruth   -0.09210    0.04189  -2.199   0.0283 * 
xtabs(~ a$toldtruth[a$DOV_PANEL=="OFF PANEL SAMPLE"] + a$demdv3[a$DOV_PANEL=="OFF PANEL SAMPLE"])
# so N is 294 not told, and 273 told truth


###
### Descriptives on undoc exposure
###

sdata <- svydesign(id=~0, weights=a$weight2, data=a)
round(svytable(~a$mstatus4, sdata),0)
#  a_cit   b_doc   c_exp d_undoc     NA
#    476     199     260     115   	 1	
#    45%     19%     25%     11%     0.1%
round(svytable(~a$fstatus4, sdata),0)
#  a_cit   b_doc   c_exp d_undoc       NA
#    452     206     251     109       32
#    43%     20%     24%     10%       3%

# so for both mothers and fathers, around 30% of those who were once undoc, are still undoc 

round(svytable(~a$mstatus4 + a$fstatus4, sdata),0)
#          a$fstatus4
#a$mstatus4 a_cit b_doc c_exp d_undoc
#   a_cit     332    54    49      20
#   b_doc      75    91    14      16
#   c_exp      30    53   168       6
#   d_undoc    16     8    19      68

# and an indicator for at least some undoc experience
a$fexp <- numeric(nrow(a))
a$fexp[a$fstatus4 =="c_exp"] <- 1
a$fexp[a$fstatus4 =="d_undoc"] <- 1
table(a$fexp)
a$mexp <- numeric(nrow(a))
a$mexp[a$mstatus4=="c_exp"] <- 1
a$mexp[a$mstatus4=="d_undoc"] <- 1
table(a$mexp)

a$someexp <- numeric(nrow(a))
a$someexp[a$fexp==1] <- 1
a$someexp[a$mexp==1] <- 1
table(a$someexp)

### code an indicator for at least some undoc experience

a$fexp <- numeric(nrow(a))
a$fexp[a$fstatus4 =="c_exp"] <- 1
a$fexp[a$fstatus4 =="d_undoc"] <- 1
table(a$fexp)
a$mexp <- numeric(nrow(a))
a$mexp[a$mstatus4=="c_exp"] <- 1
a$mexp[a$mstatus4=="d_undoc"] <- 1
table(a$mexp)

a$someexp <- numeric(nrow(a))
a$someexp[a$fexp==1] <- 1
a$someexp[a$mexp==1] <- 1
table(a$someexp)

sdata <- svydesign(id=~0, weights=a$weight2, data=a)
round(svytable(~a$someexp, sdata),0)
#  0   1 
#576 474  so that is 45%


### Knowing ppl who have been deported

round(svytable(~a$closefam_dep, sdata),0)
#   0    1 
#1002   48   so that is 5%

round(svytable(~a$distfam_dep, sdata),0)
#  0   1 
#880 170    so that is 16%

round(svytable(~a$friend_dep, sdata),0)
#  0   1 
#937 113    so that is 11%

# combine the family measures
a$close_dep <- numeric(nrow(a))
a$close_dep[a$closefam_dep==1] <- 1
a$close_dep[a$distfam_dep==1] <- 1
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
round(svytable(~a$close_dep, sdata),0)
#  0   1 
#837 213   so that is 20%

a$close_dep[a$friend_dep==1] <- 1
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
round(svytable(~a$close_dep, sdata),0)
#  0   1 
#762 288   so that is 27%



#########################################################################################
# Who knows about deportations under Obama?
#########################################################################################


## "Is obama deporting more/same/less than Bush?"
round(svytable(~a$deportq, sdata),0)
#      dk lessbush morebush  refused samebush 
#     407      165      270       11      196 
#     39%      16%      26%               19%


## Make figure showing how knowledge of deportations under Obama varies by background

# first, whether know ppl deported
a$knowdep <- numeric(nrow(a))
a$knowdep[a$closefam_dep==1 | a$distfam_dep==1 | a$friend_dep==1 | a$neigh_dep==1 | a$oth_dep==1]<-1
table(a$knowdep)

sd1 <- svydesign(id=~0, weights=a$weight2[a$knowdep==0], data=a[a$knowdep==0,])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     277      110      170        9      141 
#     40%      16%      24%               20%

sd1 <- svydesign(id=~0, weights=a$weight2[a$knowdep==1], data=a[a$knowdep==1,])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     130       55      100        2       55 
#     38%       16%     29%                16%

knowdeps <- c(30, 16, 16, 38)
knownodeps <- c(24, 20, 16, 40)

# measure of correct response: include dk as not correct, but set "refused" to NA
a$correct <- numeric(nrow(a))
a$correct[a$deportq=="morebush"] <- 1
a$correct[a$deportq=="refused"] <- NA

sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(correct ~ knowdep, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.24398    0.02717   8.981   <2e-16 ***
#knowdep      0.04979    0.04883   1.020    0.308    



# next, whether parent(s) with undoc experience
sd1 <- svydesign(id=~0, weights=a$weight2[a$someexp==0 & is.na(a$someexp)==F], data=a[a$someexp==0 & is.na(a$someexp)==F,])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     238      100      116        4      119 
#     42%      17%      20%               21%

sd1 <- svydesign(id=~0, weights=a$weight2[a$someexp==1 & is.na(a$someexp)==F], data=a[a$someexp==1 & is.na(a$someexp)==F,])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     169       65      154        8       78 
#     36%      14%      33%                17%

sglm <- svyglm(correct ~ someexp, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.20260    0.02761   7.337 4.39e-13 ***
#someexp      0.12853    0.04582   2.805  0.00512 ** 


pundoc <- c(33, 17, 14, 36)
nopundoc <- c(20, 21, 17, 42)


# next, whether spanish language media
sd1 <- svydesign(id=~0, weights=a$weight2[a$medialang=="mostspan" | a$medialang=="span" | a$medialang=="equal"], data=a[a$medialang=="mostspan" | a$medialang=="span" | a$medialang=="equal",])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     124       38       94        4       75
#     37%      11%      28%               23%

sd1 <- svydesign(id=~0, weights=a$weight2[a$medialang=="eng" | a$medialang=="mosteng"], data=a[a$medialang=="eng" | a$medialang=="mosteng",])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     283      127      176        7      122 
#     40%      18%      25%               17%

span <- c(28, 23, 11, 38)
nospan <- c(25, 17, 18, 40)

a$spanmed <- numeric(nrow(a))
a$spanmed[a$medialang=="mostspan" | a$medialang=="span" | a$medialang=="equal"] <- 1
table(a$spanmed)
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(correct ~ spanmed, design=sdata); summary(sglm)
#(Intercept)  0.24913    0.02712   9.185   <2e-16 ***
#spanmed      0.03506    0.04857   0.722    0.471  



# and political interest...
sd1 <- svydesign(id=~0, weights=a$weight2[a$polinterest=="none" | a$polinterest=="slightly"], data=a[(a$polinterest=="none" | a$polinterest=="slightly"),])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     255       62       65        9       77
#     56%      14%      14%               17%

sd1 <- svydesign(id=~0, weights=a$weight2[a$polinterest=="somewhat" | a$polinterest=="very"], data=a[(a$polinterest=="somewhat" | a$polinterest=="very"),])
t<- round(svytable(~ deportq, sd1), 0); t
#      dk lessbush morebush  refused samebush 
#     152      103      205        2      120 
#     26%      18%      35%               21%

a$polinterest_bin <- numeric(nrow(a))
a$polinterest_bin[a$polinterest=="somewhat" | a$polinterest=="very"] <- 1
table(a$polinterest_bin)
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(correct ~ polinterest_bin, design=sdata); summary(sglm)
#(Intercept)      0.14176    0.02565   5.528 4.11e-08 ***
#polinterest_bin  0.21223    0.04139   5.128 3.50e-07 ***


inter <- c(35, 21, 18, 26)
nointer <- c(14, 17, 13, 56)

depdat <- cbind(nointer, inter, nospan, span, nopundoc, pundoc, knownodeps, knowdeps)



postscript(file = "figure2.eps", horizontal=F, onefile=F, paper="special", width=6.5, height=4.5, pointsize=10)

par(mar=c(3.6, 13.9, 0.6, 6.9))
par(mgp=c(2, 0.75, 0))
barplot(depdat, horiz=T, names.arg=c(rep(NA, 8)), space=c(0.2, 0.2, 1.2, 0.2, 1.2, 0.2, 1.2, 0.2), col=c(grey(0.11), grey(0.33), grey(0.55), grey(0.77)), xlim=c(-0.5,100.5))
rect(xleft=0, ybottom=-0.5, xright=100, ytop=13, lwd=1.25)
mtext(side=2, at=1.9, line=0.25, text="Somewhat/very interested in politics", las=1)
mtext(side=2, at=0.7, line=0.25, text="Slightly/not interested in politics", las=1)
mtext(side=2, at=4.1, line=0.25, text="English-language media", las=1)
mtext(side=2, at=5.3, line=0.25, text="Spanish-language media", las=1)
mtext(side=2, at=8.7, line=0.25, text="Parent(s) with undoc. experience", las=1)
mtext(side=2, at=7.5, line=0.25, text="Parents legal migrants", las=1)
mtext(side=2, at=12.1, line=0.25, text="Acquainted with deportees", las=1)
mtext(side=2, at=10.9, line=0.25, text="Doesn't know deportees", las=1)
title(xlab="Percent of respondents")

text(119, 8.3, labels="Obama deports...", xpd=T)
legend(100, 8, legend=c('"More"', '"Same"', '"Less"', "DK"), xpd=T, fill=c(grey(0.11), grey(0.33), grey(0.55), grey(0.77)), bty="n")
text(116.5, 4.9, labels="... than Bush.", xpd=T)

dev.off()



#############################################################
# ok now experimental results, by sub-set of respondents
#############################################################

## knowledge

lm <- lm(demdv3 ~ toldtruth, data=a[a$deportq=="dk",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.61321    0.04737  12.945  < 2e-16 ***
#toldtruth   -0.23938    0.06684  -3.582 0.000424 ***
xtabs(~ a$toldtruth[a$deportq=="dk"] + a$demdv3[a$deportq=="dk"])
# so N is 106 not told truth, 107 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$deportq=="lessbush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.70787    0.05083  13.926   <2e-16 ***
#toldtruth   -0.16764    0.07230  -2.319   0.0216 *  
xtabs(~ a$toldtruth[a$deportq=="lessbush"] + a$demdv3[a$deportq=="lessbush"])
# so N is 89 not told truth, and 87 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$deportq=="samebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.49462    0.05200   9.513   <2e-16 ***
#toldtruth   -0.05018    0.07621  -0.658    0.511
xtabs(~ a$toldtruth[a$deportq=="samebush"] + a$demdv3[a$deportq=="samebush"])
# so N is 93 not told truth, and 81 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$deportq=="morebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.42857    0.04577   9.363   <2e-16 ***
#toldtruth    0.05077    0.06447   0.788    0.432  
xtabs(~ a$toldtruth[a$deportq=="morebush"] + a$demdv3[a$deportq=="morebush"])
# so N is 119 not told truth, and 121 told truth

# so those results are exactly the same as in table A3.2 and in figure 2
# this makes sense, i didn't impute these variables


## partisanship

lm <- lm(demdv3 ~ toldtruth, data=a[a$partyid=="strong_dem" | a$partyid=="notstrong_dem",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.66351    0.03340  19.864   <2e-16 ***
#toldtruth   -0.09620    0.04741  -2.029   0.0431 *  
xtabs(~ a$toldtruth[a$partyid=="strong_dem" | a$partyid=="notstrong_dem"] + a$demdv3[a$partyid=="strong_dem" | a$partyid=="notstrong_dem"])
# (that is similar to what we had before... only marginally signif)
# so N is 211 not told truth, and 208 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$partyid=="lean_dem" | a$partyid=="indep" | a$partyid=="lean_rep",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.39610    0.03878  10.215   <2e-16 ***
#toldtruth   -0.06482    0.05408  -1.199    0.232   
xtabs(~ a$toldtruth[a$partyid=="lean_dem" | a$partyid=="indep" | a$partyid=="lean_rep"] + a$demdv3[a$partyid=="lean_dem" | a$partyid=="indep" | a$partyid=="lean_rep"])
# N is 154 not told truth, 163 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$partyid=="notstrong_rep" | a$partyid=="strong_rep",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.58140    0.07474   7.779 5.06e-11 ***
#toldtruth   -0.25997    0.11901  -2.184   0.0323 *  
xtabs(~ a$toldtruth[a$partyid=="notstrong_rep" | a$partyid=="strong_rep"] + a$demdv3[a$partyid=="notstrong_rep" | a$partyid=="strong_rep"])
# N is 43 not told truth, and 28 told truth


postscript(file = "figure3.eps", horizontal=F, onefile=F, paper="special", width=5, height=3, pointsize=8.5)

par(mar=c(3.3, 15.8, 0.2, 0.2))
par(mgp=c(2, 0.75, 0))
plot(-1,1, type="n", ann=F, axes=F, xlim=c(-0.525, 0.225), ylim=c(0.58,0.97))
abline(v=0, lwd=2, col="darkgrey")
abline(v=-0.5, lwd=0.5, col="darkgrey")
abline(v=-0.4, lwd=0.5, col="darkgrey")
abline(v=-0.3, lwd=0.5, col="darkgrey")
abline(v=-0.2, lwd=0.5, col="darkgrey")
abline(v=-0.1, lwd=0.5, col="darkgrey")
abline(v=0, lwd=0.5, col="darkgrey")
abline(v=0.1, lwd=0.5, col="darkgrey")
abline(v=0.2, lwd=0.5, col="darkgrey")
box()
axis(1, at=c(-0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5), labels=c("-50", "-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50"))
title(xlab='Change in % saying Democrats "welcoming" to Latinos')

segments(-0.23938-(2*0.06684), 0.95, -0.23938+(2*0.06684), 0.95, lwd=1.5)
points(-0.23938, 0.95, pch=16, cex=1.25)
mtext(side=2, at=0.95, line=0.25, text="Don't know (n=213)", las=1)

segments(-0.16764-(2* 0.07230), 0.9, -0.16764 +(2* 0.07230), 0.9, lwd=1.5)
points(-0.16764, 0.9, pch=16, cex=1.25)
mtext(side=2, at=0.9, line=0.25, text="Obama deports less than Bush (n=176)", las=1)

segments(-0.05018-(2* 0.07621), 0.85, -0.05018 +(2* 0.07621), 0.85, lwd=1.5)
points(-0.05018, 0.85, pch=16, cex=1.25)
mtext(side=2,at=0.85, line=0.25, text="Obama deports as many as Bush (n=174)", las=1)

segments(0.05077-(2* 0.06447), 0.8, 0.05077 +(2* 0.06447), 0.8, lwd=1.5)
points(0.05077, 0.8, pch=16, cex=1.25)
mtext(side=2, at=0.8, line=0.25, text="Obama deports more than Bush (n=240)", las=1)

segments(-0.09620-(1.96* 0.04741), 0.7, -0.09620 +(1.96* 0.04741), 0.7, lwd=1.5)
points(-0.09620, 0.7, pch=15, cex=1.25)
mtext(side=2, at=0.7, line=0.25, text="Democrat (n=419)", las=1)

segments(-0.06482-(1.96* 0.05408), 0.65, -0.06482 +(1.96* 0.05408), 0.65, lwd=1.5)
points(-0.06482, 0.65, pch=15, cex=1.25)
mtext(side=2, at=0.65, line=0.25, text="Independent / none (n=317)", las=1)

segments(-0.25997-(1.96* 0.11901), 0.6, -0.25997 +(1.96* 0.11901), 0.6, lwd=1.5)
points(-0.25997, 0.6, pch=15, cex=1.25)
mtext(side=2, at=0.6, line=0.25, text="Republican (n=71)", las=1)

dev.off() 


## Ok now knowledge AND partisanship

# Democrats

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="dk",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.65385    0.06812   9.598 8.25e-16 ***
#toldtruth   -0.10283    0.09781  -1.051    0.296 
# (v similar to what i had before)
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="dk"] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="dk"])
# so N is 52 ppl not told truth, and 49 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="lessbush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.85714    0.05756  14.890  < 2e-16 ***
#toldtruth   -0.30812    0.08338  -3.695 0.000351 ***
# (so that's exactly what i had before)
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="lessbush"] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="lessbush"])
# so that is 56 people not told truth, 51 people told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="samebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.62264    0.06702   9.290 8.57e-15 ***
#toldtruth    0.01838    0.10294   0.179    0.859  
# (v similar to what i had before)
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="samebush"] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="samebush"])
# so that is 53 ppl not told truth, 39 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="morebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.50000    0.07103   7.039 1.46e-10 ***
#toldtruth    0.05882    0.09357   0.629    0.531   
# (again, v similar to what i had before)
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="morebush"] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$deportq=="morebush"])
# so that is 50 people not told truth, 68 ppl told truth


# Independents

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="dk",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.52273    0.06907   7.568 3.14e-11 ***
#toldtruth   -0.31439    0.09563  -3.288  0.00144 ** 
# (similar)
xtabs(~ a$toldtruth[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="dk"] + a$demdv3[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="dk"])
# so that is 44 ppl not told truth, 48 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="lessbush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.4074     0.0970   4.200 9.14e-05 ***
#toldtruth     0.1220     0.1299   0.939    0.352   
# (same as i had before)
xtabs(~ a$toldtruth[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="lessbush"] + a$demdv3[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="lessbush"])
# so that is 27 ppl not told truth, 34 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="samebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.29032    0.07918   3.667 0.000517 ***
#toldtruth   -0.07157    0.11110  -0.644 0.521851  
# (similar to what i had)
xtabs(~ a$toldtruth[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="samebush"] + a$demdv3[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="samebush"])
# so that is 31 ppl not told truth, 32 ppl told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="morebush",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.33333    0.06802   4.901 3.88e-06 ***
#toldtruth    0.07092    0.09822   0.722    0.472    
# (similar)
xtabs(~ a$toldtruth[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="morebush"] + a$demdv3[(a$partyid=="lean_rep" | a$partyid=="indep" | a$partyid=="lean_dem") & a$deportq=="morebush"])
# so 51 ppl not told truth, 47 ppl told truth

# (Not enough Republicans to break it down this way)



postscript(file = "figure4.eps", horizontal=F, onefile=F, paper="special", width=5, height=3, pointsize=8.5)

par(mar=c(3.3, 17.2, 0.2, 0.2))
par(mgp=c(2, 0.75, 0))
plot(-1,1, type="n", ann=F, axes=F, xlim=c(-0.5, 0.5), ylim=c(1.48,1.92))
abline(v=0, lwd=2, col="darkgrey")
abline(v=-0.5, lwd=0.5, col="darkgrey")
abline(v=-0.4, lwd=0.5, col="darkgrey")
abline(v=-0.3, lwd=0.5, col="darkgrey")
abline(v=-0.2, lwd=0.5, col="darkgrey")
abline(v=-0.1, lwd=0.5, col="darkgrey")
abline(v=0, lwd=0.5, col="darkgrey")
abline(v=0.1, lwd=0.5, col="darkgrey")
abline(v=0.2, lwd=0.5, col="darkgrey")
abline(v=0.3, lwd=0.5, col="darkgrey")
abline(v=0.4, lwd=0.5, col="darkgrey")
abline(v=0.5, lwd=0.5, col="darkgrey")
box()
axis(1, at=c(-0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5),labels=c("-50", "-40", "-30", "-20", "-10", "0", "10", "20", "30", "40", "50"))
title(xlab='Change in % saying Democrats "welcoming" to Latinos', cex.lab=0.94)

segments(-0.10283-(1.96* 0.09781), 1.9, -0.10283 +(1.96* 0.09781), 1.9, lwd=1.5)
points(-0.10283, 1.9, pch=17, cex=1.25)
mtext(side=2, at=1.9, line=0.25, text="Democrat: Don't know (n=101)", las=1)

segments(-0.30812-(1.96* 0.08338), 1.85, -0.30812 +(1.96* 0.08338), 1.85, lwd=1.5)
points(-0.30812, 1.85, pch=17, cex=1.25)
mtext(side=2, at=1.85, line=0.25, text="Democrat: Obama deports fewer (n=107)", las=1)

segments(0.01838-(1.96* 0.10294), 1.8, 0.01838 +(1.96* 0.10294), 1.8, lwd=1.5)
points(0.01838, 1.8, pch=17, cex=1.25)
mtext(side=2, at=1.8, line=0.25, text="Democrat: Obama deports as many (n=92)", las=1)

segments(0.05882-(1.96* 0.09357), 1.75, 0.05882 +(1.96* 0.09357), 1.75, lwd=1.5)
points(0.05882, 1.75, pch=17, cex=1.25)
mtext(side=2, at=1.75, line=0.25, text="Democrat: Obama deports more (n=118)", las=1)

segments(-0.31439-(1.96* 0.09563), 1.65, -0.31439 +(1.96* 0.09563), 1.65, lwd=1.5)
points(-0.31439, 1.65, pch=18, cex=1.25)
mtext(side=2, at=1.65, line=0.25, text="Independent: Don't know (n=92)", las=1)

segments(0.1220-(1.96* 0.1299), 1.6, 0.1220 +(1.96* 0.1299), 1.6, lwd=1.5)
points(0.1220, 1.6, pch=18, cex=1.25)
mtext(side=2, at=1.6, line=0.25, text="Independent: Obama deports fewer (n=61)", las=1)

segments(-0.07157-(1.96* 0.11110), 1.55, -0.07157 +(1.96* 0.11110), 1.55, lwd=1.5)
points(-0.07157, 1.55, pch=18, cex=1.25)
mtext(side=2, at=1.55, line=0.25, text="Independent: Obama deports as many (n=63)", las=1)

segments(0.07092-(1.96* 0.09822), 1.5, 0.07092 +(1.96* 0.09822), 1.5, lwd=1.5)
points(0.07092, 1.5, pch=18, cex=1.25)
mtext(side=2, at=1.5, line=0.25, text="Independent: Obama deports more (n=98)", las=1)

dev.off()


#############################################################################################
# other covariates
#############################################################################################

# political interest (among Dems)

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==0,]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.63636    0.06078   10.47   <2e-16 ***
#toldtruth   -0.17260    0.08502   -2.03   0.0444 * 
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==0] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==0])
# n = 66 not told truth, 69 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==1,]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.67586    0.03973  17.010   <2e-16 ***
#toldtruth   -0.05716    0.05680  -1.006    0.315    
xtabs(~ a$toldtruth[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==1] + a$demdv3[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem") & a$polinterest_bin==1])
# n = 145 not told truth, 139 told truth


# weaker and stronger Dems

lm <- lm(demdv3 ~ toldtruth, data=a[a$partyid=="notstrong_dem",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.60000    0.04634   12.95   <2e-16 ***
#toldtruth   -0.12419    0.06433   -1.93   0.0547 .  
# (i am going to call that p=0.05)
xtabs(~ a$toldtruth[a$partyid=="notstrong_dem"] + a$demdv3[a$partyid=="notstrong_dem"])
# n is 115 not told truth, 124 told truth

lm <- lm(demdv3 ~ toldtruth, data=a[a$partyid=="strong_dem",]); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.73958    0.04593  16.102   <2e-16 ***
#toldtruth   -0.03720    0.06724  -0.553    0.581    
xtabs(~ a$toldtruth[a$partyid=="strong_dem"] + a$demdv3[a$partyid=="strong_dem"])
# n is 96 not told truth, 84 told truth




##################################################################################
# Demographics of opt-in and panel samples
##################################################################################

# look at both weighted and unweighted results

# Age
sglm <- svyglm(PPAGE ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  23.5335     0.4007  58.730   <2e-16 ***
#DOV_PANEL2   -0.1667     0.4464  -0.373    0.709    
lm <- lm(a$PPAGE ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   24.3773     0.2118 115.097  < 2e-16 ***
#a$DOV_PANEL2  -1.2337     0.2551  -4.837 1.52e-06 ***
# so a small age difference, but only without weights; off-panel ppl a bit younger


# Parent(s) with undoc experience
sglm <- svyglm(someexp ~ DOV_PANEL, design=sdata); summary(sglm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.455275   0.051046   8.919   <2e-16 ***
#DOV_PANEL2  -0.006969   0.056589  -0.123    0.902  
lm <- lm(a$someexp ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.42638    0.02751  15.500   <2e-16 ***
#a$DOV_PANEL2  0.01837    0.03313   0.555    0.579  
# so no difference on (crucial) parental status, with or without weights


# Studying
a$anystudy <- numeric(nrow(a))
a$anystudy[a$currentstudent!="not studying"] <- 1
table(a$anystudy)
a$fullstudy <- numeric(nrow(a))
a$fullstudy[a$currentstudent=="full time"] <- 1
table(a$fullstudy)
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(anystudy ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.45963    0.05113   8.989   <2e-16 ***
#DOV_PANEL2   0.02770    0.05662   0.489    0.625   
sglm <- svyglm(fullstudy ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 0.321876   0.048038   6.700 3.39e-11 ***
#DOV_PANEL2  0.005091   0.052935   0.096    0.923    
lm <- lm(a$anystudy ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.45092    0.02767  16.298   <2e-16 ***
#a$DOV_PANEL2  0.06565    0.03332   1.971    0.049 * 
lm <- lm(a$fullstudy ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.31595    0.02634  11.994   <2e-16 ***
#a$DOV_PANEL2  0.04317    0.03172   1.361    0.174  
# so overall little evidence for edu differences


# Share with college degree
a$colleged <- numeric(nrow(a))
a$colleged[a$PPEDUC=="Bachelors degree"] <- 1
a$colleged[a$PPEDUC=="Masters degree"] <- 1
a$colleged[a$PPEDUC=="Professional or Doctorate degree"] <- 1
table(a$colleged)
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(colleged ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.08080    0.02044   3.954  8.2e-05 ***
#DOV_PANEL2   0.03430    0.02456   1.397    0.163   
lm <- lm(a$colleged ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.18712    0.02040   9.173   <2e-16 ***
#a$DOV_PANEL2 -0.03656    0.02457  -1.488    0.137    
# again, little on edu


# Parent(s) from Mexico
a$fmex <- numeric(nrow(a))
a$fmex[a$forig=="Mexico"] <- 1
table(a$fmex)
a$mmex <- numeric(nrow(a))
a$mmex[a$morig=="Mexico"] <- 1
table(a$mmex)
a$parentmex <- numeric(nrow(a))
a$parentmex[a$morig=="Mexico"] <- 1
a$parentmex[a$forig=="Mexico"] <- 1
table(a$parentmex)
sdata <- svydesign(id=~0, weights=a$weight2, data=a)
sglm <- svyglm(parentmex ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.69231    0.04554  15.201   <2e-16 ***
#DOV_PANEL2   0.01112    0.05036   0.221    0.825    
lm <- lm(a$parentmex ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.61656    0.02589  23.813  < 2e-16 ***
#a$DOV_PANEL2  0.08371    0.03118   2.685  0.00737 ** 
lm <- lm(a$fmex ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.57362    0.02648  21.665  < 2e-16 ***
#a$DOV_PANEL2  0.09765    0.03188   3.063  0.00225 ** 
lm <- lm(a$mmex ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.58589    0.02644  22.158  < 2e-16 ***
#a$DOV_PANEL2  0.08538    0.03184   2.681  0.00745 ** 
# so, without weights, this isn't balanced; off-panel ppl more mexican
# a little weird that the parents would be more mexican, but NOT more undoc
xtabs(~ a$parentmex + a$someexp)
#           a$someexp
#a$parentmex   0   1
#          0 258  84   so 25% undoc exp
#          1 331 377   so 53% undoc exp
# so those are fairly closely related... but of course not perfect

# Gender
sglm <- svyglm(female ~ DOV_PANEL, design=sdata); summary(sglm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.44217    0.05006   8.833   <2e-16 ***
#DOV_PANEL2   0.07765    0.05572   1.394    0.164  
lm <- lm(a$female ~ a$DOV_PANEL); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   0.62577    0.02702  23.162   <2e-16 ***
#a$DOV_PANEL2 -0.02079    0.03254  -0.639    0.523  
# hmm i am surprised neither of these are signif
# i recall that before imputation, there was an imbalance
# oh well, don't think i care about gender anyhow



##########################################################################################
# Check covariate balance of treated/controls
##########################################################################################

# (before weighting)

table(a$toldtruth)
#  0   1 
#534 516 

# age
lm <- lm(a$PPAGE ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  23.5805     0.1673 140.944   <2e-16 ***
#a$toldtruth  -0.1096     0.2387  -0.459    0.646   

# parents w undoc exp
lm <- lm(a$someexp ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.43071    0.02149  20.040   <2e-16 ***
#a$toldtruth  0.01696    0.03066   0.553     0.58   

lm <- lm(a$fullstudy ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.35393    0.02060  17.183   <2e-16 ***
#a$toldtruth -0.01672    0.02938  -0.569    0.569  

lm <- lm(a$colleged ~ a$toldtruth); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.162921   0.015956  10.211   <2e-16 ***
#a$toldtruth -0.002069   0.022761  -0.091    0.928    

lm <- lm(a$parentmex ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.66292    0.02029  32.667   <2e-16 ***
#a$toldtruth  0.02313    0.02895   0.799    0.425   

lm <- lm(a$female ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.60300    0.02111   28.57   <2e-16 ***
#a$toldtruth  0.01716    0.03011    0.57    0.569   

# and... partisanship, i guess
a$democrat <- numeric(nrow(a))
a$democrat[(a$partyid=="notstrong_dem" | a$partyid=="strong_dem")] <- 1
table(a$democrat)
lm <- lm(a$democrat ~ a$toldtruth); summary(lm)
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.458801   0.021578   21.26   <2e-16 ***
#a$toldtruth -0.003375   0.030781   -0.11    0.913   

lm <- lm(a$polinterest_bin ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.58801    0.02131  27.596   <2e-16 ***
#a$toldtruth  0.00307    0.03039   0.101     0.92   

lm <- lm(a$correct ~ a$toldtruth); summary(lm)
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  0.25189    0.01917  13.138   <2e-16 ***
#a$toldtruth  0.02262    0.02735   0.827    0.409  
# (so treatment status not related to whether they got the q about deportations right)

# overall, nothing here to worry about



############ END